###################################################################
# March 01 2018
# Replication File
# Rebecca Cordell
# Security-Civil Liberties Trade-offs: International Cooperation in Extraordinary Rendition
# International Interactions
# http://dvn.iq.harvard.edu/dvn/dv/internationalinteractions
###################################################################

# Load Data 
rm(list=ls())
setwd("/Users/recordel/Dropbox (ASU)/My Research/II R&R/Revisions/Final Submission/Replication Files")
data<-read.csv("Cordell-SecurityCivilLiberties-Data.csv")

# Required Packages
# install.packages("pscl")
# install.packages("lmtest")
# install.packages("MatchIt")
# install.packages("rcountrycode")
# install.packages("rworldmap")
# install.packages("dotCall64")
# install.packages("RColorBrewer")
# install.packages("cshapes")
# install.packages("maptools")
# install.packages("rgeos")
# install.packages("plyr")
# install.packages("geosphere")
# install.packages("sp")
# install.packages("ggplot2)
require("pscl")
require("lmtest")
require("MatchIt")
require("countrycode")
require("rworldmap")  
require("dotCall64")
require("RColorBrewer")
require("cshapes")
require("maptools")
require("rgeos")
require("plyr")
require("geosphere")
require("sp")
require("ggplot2")
options(scipen=999)

###################################################################
# Table 1: Descriptive Statistics of Independent and Control 
# Variables
###################################################################
# N
hr_sim_9100 <- data[!is.na(data$hr_sim_9100),]
all_sim_9100 <- data[!is.na(data$all_sim_9100),]
centrendcorr_log <- data[!is.na(data$centrendcorr),]
alliance <- data[!is.na(data$alliance),]
regime <- data[!is.na(data$regime),]
terr_9100_log <- data[!is.na(data$terr_9100_log),]
ustrade_log <- data[!is.na(data$ustrade_log),]
usaid_log <- data[!is.na(data$usaid_log),]
pop_log <- data[!is.na(data$pop_log),]
partyideo <- data[!is.na(data$partyideo),]
n_hrsim9100<-length(hr_sim_9100$hr_sim_9100)
n_allsim9100<-length(all_sim_9100$all_sim_9100)
n_centrendcorrlog<-length(centrendcorr_log$centrendcorr_log)
n_alliance<-length(alliance$alliance)
n_regime<-length(regime$regime)
n_terr9100log<-length(terr_9100_log$terr_9100_log)
n_ustradelog<-length(ustrade_log$ustrade_log)
n_usaidlog<-length(usaid_log$usaid_log)
n_poplog<-length(pop_log$pop_log)
n_partyideo<-length(partyideo$partyideo)

# Standard Deviation
sd_hrsim9100<-sd(hr_sim_9100$hr_sim_9100) 
sd_allsim9100<-sd(all_sim_9100$all_sim_9100)
sd_centrendcorrlog<-sd(centrendcorr_log$centrendcorr_log)
sd_alliance<-sd(alliance$alliance)
sd_regime<-sd(regime$regime)
sd_terr9100log<-sd(terr_9100_log$terr_9100_log)
sd_ustradelog<-sd(ustrade_log$ustrade_log)
sd_usaidlog<-sd(usaid_log$usaid_log)
sd_poplog<-sd(pop_log$pop_log) 
sd_partyideo<-sd(partyideo$partyideo)

# Mean
mean_hrsim9100<-mean(hr_sim_9100$hr_sim_9100) 
mean_allsim9100<-mean(all_sim_9100$all_sim_9100)
mean_centrendcorrlog<-mean(centrendcorr_log$centrendcorr_log)
mean_alliance<-mean(alliance$alliance)
mean_regime<-mean(regime$regime)
mean_terr9100log<-mean(terr_9100_log$terr_9100_log)
mean_ustradelog<-mean(ustrade_log$ustrade_log)
mean_usaidlog<-mean(usaid_log$usaid_log)
mean_poplog<-mean(pop_log$pop_log) 
mean_partyideo<-mean(partyideo$partyideo)

# Min
min_hrsim9100<-min(hr_sim_9100$hr_sim_9100) 
min_allsim9100<-min(all_sim_9100$all_sim_9100)
min_centrendcorrlog<-min(centrendcorr_log$centrendcorr_log)
min_alliance<-min(alliance$alliance)
min_regime<-min(regime$regime)
min_terr9100log<-min(terr_9100_log$terr_9100_log)
min_ustradelog<-min(ustrade_log$ustrade_log)
min_usaidlog<-min(usaid_log$usaid_log)
min_poplog<-min(pop_log$pop_log) 
min_partyideo<-min(partyideo$partyideo)

# Max
max_hrsim9100<-max(hr_sim_9100$hr_sim_9100) 
max_allsim9100<-max(all_sim_9100$all_sim_9100)
max_centrendcorrlog<-max(centrendcorr_log$centrendcorr_log)
max_alliance<-max(alliance$alliance)
max_regime<-max(regime$regime)
max_terr9100log<-max(terr_9100_log$terr_9100_log)
max_ustradelog<-max(ustrade_log$ustrade_log)
max_usaidlog<-max(usaid_log$usaid_log)
max_poplog<-max(pop_log$pop_log) 
max_partyideo<-max(partyideo$partyideo)

# Create table
sumstats_hrsim9100<-cbind(n_hrsim9100, mean_hrsim9100, sd_hrsim9100, min_hrsim9100, max_hrsim9100)
sumstats_allsim9100<-cbind(n_allsim9100, mean_allsim9100, sd_allsim9100, min_allsim9100, max_allsim9100)
sumstats_centrendcorrlog<-cbind(n_centrendcorrlog, mean_centrendcorrlog, sd_centrendcorrlog, min_centrendcorrlog, max_centrendcorrlog)
sumstats_alliance<-cbind(n_alliance, mean_alliance, sd_alliance, min_alliance, max_alliance)
sumstats_regime<-cbind(n_regime, mean_regime, sd_regime, min_regime, max_regime)
sumstats_terr9100log<-cbind(n_terr9100log, mean_terr9100log, sd_terr9100log, min_terr9100log, max_terr9100log)
sumstats_ustradelog<-cbind(n_ustradelog, mean_ustradelog, sd_ustradelog, min_ustradelog, max_ustradelog)
sumstats_usaidlog<-cbind(n_usaidlog, mean_usaidlog, sd_usaidlog, min_usaidlog, max_usaidlog)
sumstats_poplog<-cbind(n_poplog, mean_poplog, sd_poplog, min_poplog, max_poplog)
sumstats_partyideo<-cbind(n_partyideo, mean_partyideo, sd_partyideo, min_partyideo, max_partyideo)
table1<-rbind(sumstats_hrsim9100, sumstats_allsim9100, sumstats_centrendcorrlog, sumstats_alliance, sumstats_regime, sumstats_terr9100log, sumstats_ustradelog, sumstats_usaidlog, sumstats_poplog, sumstats_partyideo)
table1<-as.data.frame(table1)
names(table1) <- gsub("_hrsim9100", "", names(table1))
table1<-round(table1, digits=2)
table1

###################################################################
# Table 2: Probit Regression, Cooperation in RDI
###################################################################
# Model 1 Baseline Model Human Rights Dissimilarity
table2_model1<- glm(osf_coop ~ hr_sim_9100, family=binomial(link="probit"), data=data)
# Estimate, Std Error, AIC and N (N = 2 observations deleted due to missingness 169-2=197)
print(summary(table2_model1, robust=TRUE), digits=2)
# Log Likelihood
round(logLik(table2_model1), digits=2)
# Pseudo R2: McFadden
round(pR2(table2_model1), digits=2)
# LR chi2 and Prob>chi2 
round(lrtest(table2_model1), digits=2)

# Model 2 Full Model Human Rights Dissimilarity
table2_model2<- glm(osf_coop ~ hr_sim_9100 + all_sim_9100 + centrendcorr_log + alliance + regime + terr_9100_log + ustrade_log + usaid_log + pop_log + partyideo, family=binomial(link="probit"), data=data)
# Estimate, Std Error, AIC and N (N = 14 observations deleted due to missingness 169-14=155)
print(summary(table2_model2, robust=TRUE),digits=1)
# Log Likelihood
round(logLik(table2_model2), digits=2)
# Pseudo R2: McFadden
round(pR2(table2_model2), digits=2)
# LR chi2 and Prob>chi2 
round(lrtest(table2_model2), digits=2)

###################################################################
# Table 3: Estimating Model Accuracy with Probit Regression, 
# Cooperation in RDI – Using Matched Sample
###################################################################
# Create matched sample using propensity score matching
# Treatment group: hr_sim_9100 <= 1.258081 (assigned a value of 1, and 0 otherwise)
matching<-na.omit(data)
matching<-matching[,c(2, 4:15)]
matched_model <- matchit(treatment ~ all_sim_9100 + centrendcorr_log + alliance + regime + terr_9100_log + ustrade_log + usaid_log + pop_log + partyideo,
                     method = "nearest", data = matching)
matched_data <- match.data(matched_model)

# Model 1 Baseline Model Human Rights Dissimilarity
table3_model1<- glm(osf_coop ~ hr_sim_9100, family=binomial(link="probit"), data=matched_data)
# Estimate, Std Error and AIC (N = 126 obs)
print(summary(table3_model1, robust=TRUE), digits=2)
# Log Likelihood
round(logLik(table3_model1), digits=2)
# Pseudo R2: McFadden
round(pR2(table3_model1), digits=2)
# LR chi2 and Prob>chi2 
round(lrtest(table3_model1), digits=2)

# Model 2 Full Model Human Rights Dissimilarity
table3_model2<- glm(osf_coop ~ hr_sim_9100 + all_sim_9100 + centrendcorr_log + alliance + regime + terr_9100_log + ustrade_log + usaid_log + pop_log + partyideo, family=binomial(link="probit"), data=matched_data)
# Estimate, Std Error and AIC (N = 126 obs)
print(summary(table3_model2, robust=TRUE),digits=1)
# Log Likelihood
round(logLik(table3_model2), digits=2)
# Pseudo R2: McFadden
round(pR2(table3_model2), digits=2)
# LR chi2 and Prob>chi2 
round(lrtest(table3_model2), digits=2)

###################################################################
# Figure 1: International Cooperation in Rendition, Secret 
# Detention and Interrogation (Alliances)
###################################################################
# Create subsample of states that cooperated
figure1<-data[data $ osf_coop== 1,]
# Assign ISO country code to map data:
figure1 $ isocode <- countrycode(figure1 $ statenme, origin = "country.name", destination = "iso3c")
# Join data and country map:
map.data <- joinCountryData2Map(        
  figure1,                         
  joinCode = "ISO3",                    
  nameJoinColumn = "isocode",           
  verbose = TRUE)
# Create colour palette
colourPalette <- palette(c("turquoise1", "blue1"))
# Plot map
mapParams <- mapCountryData(map.data,                                     
                            nameColumnToPlot="alliance",                 
                            addLegend=FALSE,  
                            colourPalette=colourPalette,
                            catMethod="categorical",                                                 
                            missingCountryCol="white", 
                            mapTitle="", 
                            borderCol="black")  

###################################################################
# Figure 2: Central Rendition Transit Corridor Control Variable
###################################################################
# Compile shape file # 
cshp(date=NA, useGW=TRUE)
cshp.data <- cshp()
wmap <- cshp(date=as.Date("2000-09-11"))
wmap_df <- fortify(wmap)
plot(wmap,col="#f2f2f2", fill=TRUE, bg="white", lwd=0.5)
# Plot closest airport within a country to the central rendition transit corridor
points(data$long,data$lat, pch=10, cex=0.2, col="blue")
# Draw a great circle (shortest flight path from Washington Dulles International Airport, U.S. to Kabul International Airport, Afghanistan)
lat_us <- 38.9445 # Washington Dulles International Airport (KIAD)
long_us <- -77.4558029 # Washington Dulles International Airport (KIAD)
lat_af <- 34.5658989 # Kabul International Airport (OAKB)
long_af <- 69.2123032 # Kabul International Airport (OAKB)
inter <- gcIntermediate(c(long_us, lat_us), c(long_af, lat_af), n=50, addStartEnd=TRUE)
lines(inter,lwd=2, col="red")

###################################################################
# Figure 3: Predicted Probabilities with 95% Confidence Intervals 
# (Model 2)
###################################################################
# Predicted Probabilities holding all Controls at their means #
# Mean
predprob<-na.omit(data)
predprob<-predprob[,c(4:14)]
predprob$all_sim_9100<-mean(predprob$all_sim_9100)
predprob$centrendcorr_log<-mean(predprob$centrendcorr_log)
predprob$alliance<-mean(predprob$alliance)
predprob$regime<-mean(predprob$regime)
predprob$terr_9100_log<-mean(predprob$terr_9100_log)
predprob$ustrade_log<-mean(predprob$ustrade_log)
predprob$usaid_log<-mean(predprob$usaid_log)
predprob$pop_log<-mean(predprob$pop_log) 
predprob$partyideo<-mean(predprob$partyideo)

figure3_table1_model2<- glm(osf_coop ~ hr_sim_9100 + all_sim_9100 + centrendcorr_log + alliance + regime + terr_9100_log + ustrade_log + usaid_log + pop_log + partyideo, family=binomial(link="probit"), data=predprob)
predprob$pred.prob <- predict(figure3_table1_model2, newdata=predprob, type="response")

# Adjust the x axis labels #
x1 <- seq(0, 5, 0.01)
x2 <- seq(0, 5, 0.01)
x3 <- seq(0, 5, 0.01)
x4 <- seq(0, 5, 0.01)
x5 <- seq(0, 5, 0.01)
x6 <- seq(0, 5, 0.01)
x7 <- seq(0, 5, 0.01)
x8 <- seq(0, 5, 0.01)
x9 <- seq(0, 5, 0.01)
x10 <- seq(0, 5, 0.01)
hr_sim_91000 <- seq(0, 5, 0.01)
newdf <- data.frame(x = seq(0, 1, length.out = 501))
newdf$pout_probit <- predict(figure3_table1_model2, newdata=data.frame(hr_sim_9100=x1,all_sim_9100=x2, centrendcorr_log=x3, alliance=x4, regime=x5, terr_9100_log=x6, ustrade_log=x7, usaid_log=x8, pop_log=x9, partyideo=x10), se.fit = TRUE, type = "response")$fit
newdf$pse_probit <- predict(figure3_table1_model2, newdata=data.frame(hr_sim_9100=x1,all_sim_9100=x2, centrendcorr_log=x3, alliance=x4, regime=x5, terr_9100_log=x6, ustrade_log=x7, usaid_log=x8, pop_log=x9, partyideo=x10), se.fit = TRUE, type = "response")$se.fit
newdf$pupper_probit <- newdf$pout_probit + (1.96 * newdf$pse_probit)
newdf$plower_probit <- newdf$pout_probit - (1.96 * newdf$pse_probit)
plot(predprob$osf_coop ~ predprob$hr_sim_9100, col = NULL, bg = rgb(0, 0, 0, 0.5), pch = 21, las=1, ylab="Predicted Probabilities", xlab="Human Rights Dissimilarity", par(las=3))
with(newdf, lines(pout_probit ~ x1, type = "l", lwd = 2, col="darkblue" ))
with(newdf, lines(pupper_probit ~ x1, type = "l", lty = 2))
with(newdf, lines(plower_probit ~ x1, type = "l", lty = 2))

###################################################################
# Appendix 1: Categories of participation in the RDI program
###################################################################
appendix1<-data[,c(2, 18:21)]
appendix1[appendix1=="1"]<-"Yes"
appendix1[appendix1=="0"]<-"No"
appendix1

###################################################################
# Appendix 2: List of UNGA human rights resolutions from 1991-2000
###################################################################
# See Vote Description Data published here for the subjects of all UNGA human rights resolutions from 1991-2000: https://hdl.handle.net/1902.1/12379 
# Erik Voeten "Data and Analyses of Voting in the UN General Assembly" Routledge Handbook of International Organization, edited by Bob Reinalda (published May 27, 2013). Available at SSRN: http://ssrn.com/abstract=2111149

###################################################################
# Appendix 3: Matched Sample used to estimate treatment effect of 
# Human Rights Similarity on Cooperation in the RDI program
###################################################################
treatment<-matched_data[matched_data $ treatment== 1,]
control<-matched_data[matched_data $ treatment== 0,]
View(treatment)
treatment<-as.data.frame(treatment[,c(1)])
control<-as.data.frame(control[,c(1)])
appendix3<-cbind(treatment, control)
appendix3

###################################################################
# Appendix 4: Probit Regression, Cooperation in RDI – Using Economic 
# Development Dissimilarity in replacement of Human Rights Dissimilarity
###################################################################
# Model 1 Baseline Model Human Rights Dissimilarity
appendix4_model1<- glm(osf_coop ~ ed_sim_9100, family=binomial(link="probit"), data=data)
# Estimate, Std Error, AIC and N (N = 2 observations deleted due to missingness 169-2=197)
print(summary(appendix4_model1, robust=TRUE), digits=2)
# Log Likelihood
round(logLik(appendix4_model1), digits=2)
# Pseudo R2: McFadden
round(pR2(appendix4_model1), digits=2)
# LR chi2 and Prob>chi2 
round(lrtest(appendix4_model1), digits=2)

# Model 2 Full Model Human Rights Dissimilarity
appendix4_model2<- glm(osf_coop ~ ed_sim_9100 + all_sim_9100 + centrendcorr_log + alliance + regime + terr_9100_log + ustrade_log + usaid_log + pop_log + partyideo, family=binomial(link="probit"), data=data)
# Estimate, Std Error, AIC and N (N = 14 observations deleted due to missingness 169-14=155)
print(summary(appendix4_model2, robust=TRUE),digits=1)
# Log Likelihood
round(logLik(appendix4_model2), digits=2)
# Pseudo R2: McFadden
round(pR2(appendix4_model2), digits=2)
# LR chi2 and Prob>chi2 
round(lrtest(appendix4_model2), digits=2)

###################################################################
# P.5 "However, this arguments fails to fully account for international cooperation in RDI operations as only 28% of countries that participated had formal alliances with the U.S.."
###################################################################
coop_alliance<-alliance[alliance$osf_coop== 1,] # 53 observations
coop_alliance$alliance<-as.factor(coop_alliance$alliance)
summary(coop_alliance$alliance) # 15 observations
round(15/53*100, digits=0)

###################################################################
# P.5 "However, this also fails to explain international cooperation in RDI, since as many as 40% of participating states were democracies, including established democracies such as Denmark and Sweden."
###################################################################
coop_regime<-regime[regime$osf_coop== 1,] # 52 observations
median(regime$regime)
regime_median<-coop_regime[coop_regime$regime>= 0.5529295,] # 26 observations
round(26/52*100, digits=0)

###################################################################
# P.6 "However, 43% of the countries that participated in RDI operations actually faced low levels of terrorism threat, with participating states such as Finland and Portugal not experiencing a single terrorism attack in the ten years that preceded 9/11."
###################################################################
coop_terr_9100_log<-terr_9100_log[terr_9100_log$osf_coop== 1,] # 53 observations
median(terr_9100_log$terr_9100_log)
terr_9100_log_median<-coop_terr_9100_log[coop_terr_9100_log$terr_9100_log<= 1.791759,] # 21 observations
round(21/53*100, digits=0)

###################################################################
# P.7 "For example, between 40%-49% of the countries that cooperated had low levels of U.S. trade and aid, and 66% had a large population size above the global median."
###################################################################
coop_ustrade_log<-ustrade_log[ustrade_log$osf_coop== 1,] # 53 observations
median(ustrade_log$ustrade_log)
ustrade_log_median<-coop_ustrade_log[coop_ustrade_log$ustrade_log<= 8.847966,] # 29 observations
round(29/53*100, digits=0)

coop_usaid_log<-usaid_log[usaid_log$osf_coop== 1,] # 53 observations
median(usaid_log$usaid_log)
usaid_log_median<-coop_usaid_log[coop_usaid_log$usaid_log<= 6.038888,] # 33 observations
round(33/53*100, digits=0)

coop_pop_log<-pop_log[pop_log$osf_coop== 1,] # 53 observations
median(pop_log$pop_log)
pop_log_median<-coop_pop_log[coop_pop_log$pop_log>= 9.044562,] # 85 observations
round(34/53*100, digits=0)

###################################################################
# P.7 "However, we find quite the opposite, as 81% of the countries that participated in RDI operations had left-of-center governments (including countries such as Canada and the UK)."
###################################################################
coop_partyideo<-partyideo[partyideo$osf_coop== 1,] # 53 observations
coop_partyideo$partyideo<-as.factor(coop_partyideo$partyideo)
summary(coop_partyideo$partyideo) # 46 observations
round(46/53*100, digits=0)

###################################################################
# P.13 "This is the case for 31% of the observations (53) but not for the remaining 69% (116)."
###################################################################
data$osf_coop<-as.factor(data$osf_coop)
summary(data$osf_coop) # 46 observations
round(53/169*100, digits=0)
round(116/169*100, digits=0)

###################################################################
# P.14 "The index is continuous and has a scale from 0.090 to 2.194. States with lower values (greater similarity) such as Israel have more similar human rights preferences to the U.S. and states with higher values (greater dissimilarity) such as China have less similar human rights preferences to the U.S.."
###################################################################
round(min(hr_sim_9100$hr_sim_9100), digits=2) # minimum
round(max(hr_sim_9100$hr_sim_9100), digits=2) # maximum

###################################################################
# P.14 "The median value is 31.250 (for example, Senegal), the minimum value is 0.090 (Israel), the maximum value is 2.194 (India) and the standard deviation across the sample is 0.440."
###################################################################
round(median(hr_sim_9100$hr_sim_9100), digits=2) # median
round(min(hr_sim_9100$hr_sim_9100), digits=2) # minimum
round(max(hr_sim_9100$hr_sim_9100), digits=2) # maximum
round(sd(hr_sim_9100$hr_sim_9100), digits=2) # standard deviation

###################################################################
# P. 16 "The index is continuous and has a scale from 0.520 to 4.518."
###################################################################
round(min(all_sim_9100$all_sim_9100), digits=2) # minimum
round(max(all_sim_9100$all_sim_9100), digits=2) # maximum

###################################################################
# P. 26 "When the independent variable is held at its mean (for example, Cyprus and Portugal), the probability of cooperation in RDI is 43%."
###################################################################
mean(hr_sim_9100$hr_sim_9100) # mean
mean <- subset(hr_sim_9100, hr_sim_9100 <= 1.258081, 
                     select=c(ccode, osf_coop, hr_sim_9100))
mean_predprob<- glm(osf_coop ~ hr_sim_9100, family=binomial(link="probit"), data=mean)
mean_predprob$pred.prob <- predict(mean_predprob, newdata=mean, type="response")  
round(mean(mean_predprob$pred.prob), digits=2)

###################################################################
# P. 26 "When the independent variable is held at the 95th percentile (for example, the Australia and Ireland), the probability of cooperation in RDI increases to 56%."
###################################################################
quantile(hr_sim_9100$hr_sim_9100, 0.05) # 95th percentile
quant95th <- subset(hr_sim_9100, hr_sim_9100 <= 0.495425, 
                     select=c(ccode, osf_coop, hr_sim_9100))
quant95th_predprob<- glm(osf_coop ~ hr_sim_9100, family=binomial(link="probit"), data=quant95th)
quant95th_predprob$pred.prob <- predict(quant95th_predprob, newdata=quant_95th, type="response")                      
round(mean(quant95th_predprob$pred.prob), digits=2)
